!****************************************************************************************
SUBROUTINE IRFS2
! compute IRFs
! NOTE: Set na2=501, nb2=3501

USE GLOBAL

IMPLICIT NONE

INTEGER::nirf
PARAMETER (nirf=200)

REAL(8),DIMENSION(nirf,nmc)::at,pt,ht,ctt,gnt,bt,acct,deft,rt,wt,taut,bort
REAL(8),DIMENSION(nirf,nmc)::yt,ynt,ytt,cnt,ct

REAL(8)::at0,gnt0,ctt0,ht0,pt0,wt0,bt0,taut0,rt0,b00,a_old,a_oldss,a0,a0ss 
REAL(8)::at0ss,gnt0ss,ctt0ss,ht0ss,pt0ss,wt0ss,bt0ss,taut0ss,rt0ss  
REAL(8)::cnt0ss,ct0ss,ynt0ss,ytt0ss,yt0ss,yttv0ss,ytv0ss
REAL(8)::at2(nt2),gnt2(nt2),ctt2(nt2),ht2(nt2),cnt2(nt2),pt2(nt2),wt2(nt2)
REAL(8)::bt2(nt2),rt2(nt2),taut2(nt2),yt2(nt2),ytv2(nt2),ynt2(nt2),ytt2(nt2),yttv2(nt2),ct2(nt2),dummy1(nirf),bss

INTEGER::i_at(nirf,nmc),i_bt(nirf,nmc),i_at0(nirf,nmc),i_a0,bort0,acct0,acct02,deft0,bort0ss,acct0ss,deft0ss
INTEGER::i,j,k,t,i_a,i_b,i_ass,i_bss,i_aut0(1),i_at0ss(nirf,nmc),i_b0(1),i_b0ss,i_b0ss0(1),i_at_initial

REAL::eps_a_vec(nirf)
REAL(8)::eps_a

1002 FORMAT(1000F7.2)

seed=1234658753
CALL RNSET (seed)

i_a0 = 1
a0ss = mua/(1d0-rhoa) !+1d0*sigmaa  !start from unconditional mean

!SS debt indices for fine grids (nb2=3501)
!flex wages: 849
!baseline: 745
!risk-free debt: 7883

!i_bss = 745                 !steady-state debt
!bss  = 0.75d0*bgrid2(i_bss) 
!b00  = bss 
!i_b0 = LOCATE2(bgrid2,b00)

i_at_initial = LOCATE2(agrid2,a0ss)
i_at0ss(1,:) = i_at_initial

i_b0ss0=MINLOC((bgrid2(bp1(i_at_initial,:))-bgrid2)**2d0)
i_bss=i_b0ss0(1)
PRINT*,'i_at_initial,index_b',i_at_initial,i_bss

!low debt:  0.75*bss
!high debt: bss

bss  = 0.75d0*bgrid2(i_bss)         

!b00: initial debt in IRFs
b00  = bss 
i_b0 = LOCATE2(bgrid2,b00)

!a0: initial yT in IRFs
a0 = a0ss-1d0*sigmaa
i_at_initial = LOCATE2(agrid2,a0)    
i_at0(1,:) = i_at_initial
 

DO j=1,nmc

a_oldss  = a0ss
a_old    = a0 

i_b   = i_b0(1)
i_bss = i_b0(1) 

acct0 = 1   !access in initial period
deft0 = 0
bort0 = 1
    
acct0ss= 1   !access in initial period
deft0ss= 0
bort0ss= 1

CALL r8vec_uniform_01(nirf,seed,dummy1)    !reentry  
CALL RNNOR(nirf,eps_a_vec)

dummy1(1)= 1d0        !no reentry after default

IF (def1(i_at0(1,j),i_b).EQ.0) THEN
    PRINT*,'Default occurs in initial state',j,i_at0(1,j),i_b
    PAUSE
    STOP
ENDIF
    
DO i=1,nirf

    eps_a = sigmaa*DBLE(eps_a_vec(i))
    
    IF (i.GE.2) THEN
        i_at0(i,j) = LOCATE2(agrid2,rhoa*a_old+eps_a)
        a_old    = rhoa*a_old+eps_a
        i_at0ss(i,j) = LOCATE2(agrid2,rhoa*a_oldss+eps_a)
        a_oldss    = rhoa*a_oldss+eps_a
    ENDIF
    
    i_a = i_at0(i,j)
    i_ass=i_at0ss(i,j)

    acct02 = acct0
    deft0  = 0
    bort0  = 1
    deft0ss= 0
    bort0ss= 1

    IF ((acct0.EQ.0).AND.(dummy1(i)>phi)) THEN
        bort0 = 0 !bort(i)=0
        acct0 = 0 !acct(i+1)=0
    ELSEIF ((acct0.EQ.0).AND.(dummy1(i).LE.phi)) THEN
        bort0 = 1 !bort(i)=0
        acct0 = 1 !acct(i+1)=0
    ENDIF    
    
    IF (acct0.EQ.1) THEN 
        IF (def1(i_a,i_b).EQ.1d0) THEN
            bort0 = 1 !bort(i)=1
            acct0 = 1 !acct(i+1)=1
        ELSE
            bort0 = 0  !bort(i)=0
            acct0 = 0  !acct(i+1)=0
            deft0 = 1  !deft(i)=1
        ENDIF
    ENDIF
    
    IF (bort0.EQ.1) THEN    
        at0  = DEXP(agrid2(i_a))
        gnt0 = gn1(i_a,i_b)
        ctt0 = ct1(i_a,i_b)
        ht0  = h1(i_a,i_b)
        pt0  = p1(i_a,i_b)
        wt0  = w1(i_a,i_b)
        taut0= tau1(i_a,i_b)
        i_b   = bp1(i_a,i_b)
        bt0   = bgrid2(i_b)/(lambda+r)
        acct0 = 1
        rt0   = (1d0+q1(i_a,i_b)**(-1d0)-lambda)/(1d0+r)-1d0

    ELSEIF (bort0.EQ.0) THEN
        at0  = DEXP(agrid2(i_a))
        gnt0 = gnaut1(i_a)
        ctt0 = ctaut1(i_a)
        ht0  = haut1(i_a)
        pt0  = paut1(i_a)
        wt0  = waut1(i_a)
        taut0= tauaut1(i_a)
        rt0  = -1d0
        i_b  = i_bzero2
    ENDIF

    IF ((acct0ss.EQ.0).AND.(dummy1(i)>phi)) THEN
        bort0ss = 0 !bort(i)=0
        acct0ss = 0 !acct(i+1)=0
    ELSEIF ((acct0ss.EQ.0).AND.(dummy1(i).LE.phi)) THEN
        bort0ss = 1 !bort(i)=0
        acct0ss = 1 !acct(i+1)=0
    ENDIF
    
    IF (acct0ss.EQ.1) THEN 
        IF (def1(i_ass,i_bss).EQ.1d0) THEN
            bort0ss = 1 !bort(i)=1
            acct0ss = 1 !acct(i+1)=1
        ELSE
            bort0ss = 0  !bort(i)=0
            acct0ss = 0  !acct(i+1)=0
            deft0ss = 1  !deft(i)=1
        ENDIF
    ENDIF

    IF (bort0ss.EQ.1) THEN    
        at0ss  = DEXP(agrid2(i_ass))
        gnt0ss = gn1(i_ass,i_bss)
        ctt0ss = ct1(i_ass,i_bss)
        ht0ss  = h1(i_ass,i_bss)
        pt0ss  = p1(i_ass,i_bss)
        wt0ss  = w1(i_ass,i_bss)
        taut0ss= tau1(i_ass,i_bss)

        i_bss   = bp1(i_ass,i_bss)
        bt0ss   = bgrid2(i_bss)/(lambda+r)
        acct0ss = 1
        rt0ss   = (1d0+q1(i_ass,i_bss)**(-1d0)-lambda)/(1d0+r)-1d0

    ELSEIF (bort0ss.EQ.0) THEN
        at0ss  = DEXP(agrid2(i_ass))
        gnt0ss = gnaut1(i_ass)
        ctt0ss = ctaut1(i_ass)
        ht0ss  = haut1(i_ass)
        pt0ss  = paut1(i_ass)
        wt0ss  = waut1(i_ass)
        taut0ss= tauaut1(i_ass)
        rt0ss  = -1d0
        i_bss  = i_bzero2
    ENDIF

    i_at(i,j)= i_a
    at(i,j)  = at0-at0ss
    gnt(i,j) = gnt0-gnt0ss
    ctt(i,j) = ctt0-ctt0ss
    ht(i,j)  = ht0-ht0ss
    pt(i,j)  = pt0-pt0ss
    wt(i,j)  = wt0-wt0ss
    i_bt(i,j)= i_b
    bt(i,j)  = bgrid2(i_b)-bgrid2(i_bss)
    rt(i,j)  = rt0-rt0ss
    taut(i,j)= taut0-taut0ss
    deft(i,j)= deft0-deft0ss
    acct(i,j)= acct02-acct0ss
    bort(i,j)= DBLE(bort0)-DBLE(bort0ss)


    bt(i,j)  = bt(i,j)/(lambda+r)
    cnt0ss = ht0ss**alpha-gnt0ss

    IF (mu.EQ.0d0) THEN
        ct0ss = (ctt0ss**(omegac))*(cnt0ss**(1d0-omegac))
    ELSE
        ct0ss = (omegac*ctt0ss**(-mu)+(1d0-omegac)*cnt0ss**(-mu))**(-1d0/mu)
    ENDIF
    
    ynt0ss = ht0ss**alpha
    ytt0ss = at0ss
    yt0ss  = ytt0ss + pt0ss*ynt0ss 
    yttv0ss= ytt0ss         
    ytv0ss = yt0ss 
    
    cnt(i,j) = ht0**alpha-gnt0-cnt0ss

    IF (mu.EQ.0d0) THEN
        ct(i,j) = (ctt(i,j)**(omegac))*(cnt(i,j)**(1d0-omegac))-ct0ss
    ELSE
        ct(i,j) = (omegac*ctt(i,j)**(-mu)+(1d0-omegac)*cnt(i,j)**(-mu))**(-1d0/mu)-ct0ss
    ENDIF
    
    ynt(i,j) = ht0**alpha -ynt0ss
    ytt(i,j) = at0 -ytt0ss
    yt(i,j)  = ytt(i,j) + pt0*ynt(i,j)-yt0ss 

ENDDO  !i
ENDDO  !j


OPEN(1,FILE='output//irf//irfp.txt')
OPEN(2,FILE='output//irf//irfh.txt')
OPEN(3,FILE='output//irf//irfb.txt')
OPEN(4,FILE='output//irf//irfct.txt')
OPEN(5,FILE='output//irf//irfr.txt')
OPEN(7,FILE='output//irf//irfgn.txt')
OPEN(8,FILE='output//irf//irftau.txt')    
OPEN(9,FILE='output//irf//irfa.txt')      
OPEN(11,FILE='output//irf//irfcn.txt')    
OPEN(12,FILE='output//irf//irfww.txt')    
OPEN(13,FILE='output//irf//irfbor.txt') 

DO i=1,nirf
    WRITE(1,'(F15.10,F15.10,F15.10)') quantile(pt(i,:),0.5d0), quantile(pt(i,:),0.25d0), quantile(pt(i,:),0.75d0) 
    WRITE(2,'(F15.10,F15.10,F15.10)') quantile(ht(i,:),0.5d0), quantile(ht(i,:),0.25d0), quantile(ht(i,:),0.75d0)
    WRITE(3,'(F15.10,F15.10,F15.10)') quantile(bt(i,:),0.5d0), quantile(bt(i,:),0.25d0), quantile(bt(i,:),0.75d0)
    WRITE(4,'(F15.10,F15.10,F15.10)') quantile(ctt(i,:),0.5d0), quantile(ctt(i,:),0.25d0), quantile(ctt(i,:),0.75d0)
    WRITE(5,'(F15.10,F15.10,F15.10)') quantile(rt(i,:),0.5d0), quantile(rt(i,:),0.25d0), quantile(rt(i,:),0.75d0)
    WRITE(7,'(F15.10,F15.10,F15.10)') quantile(gnt(i,:),0.5d0), quantile(gnt(i,:),0.25d0), quantile(gnt(i,:),0.75d0)
    WRITE(8,'(F15.10,F15.10,F15.10)') quantile(taut(i,:),0.5d0), quantile(taut(i,:),0.25d0), quantile(taut(i,:),0.75d0) 
    WRITE(9,'(F15.10,F15.10,F15.10)') quantile(at(i,:),0.5d0), quantile(at(i,:),0.25d0), quantile(at(i,:),0.75d0)
    WRITE(11,'(F15.10,F15.10,F15.10)')quantile(cnt(i,:),0.5d0), quantile(cnt(i,:),0.25d0), quantile(cnt(i,:),0.75d0)
    WRITE(12,'(F15.10,F15.10,F15.10)') quantile(wt(i,:),0.5d0), quantile(wt(i,:),0.25d0), quantile(wt(i,:),0.75d0)
    WRITE(13,'(F15.10,F15.10,F15.10)') quantile(bort(i,:),0.5d0), quantile(bort(i,:),0.25d0), quantile(bort(i,:),0.75d0)
ENDDO

CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(11)
CLOSE(12)
CLOSE(13)
    
PRINT*, ' '
PRINT*, 'IMPUlSE RESPONSES EXPORTED'
PRINT*, ' '
    
END SUBROUTINE IRFS2
